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In previous work, we analyzed the linear and nonlinear stability of static, spherically symmetric 
wormhole solutions to Einstein's field equations coupled to a massless ghost scalar field. Our anal- 
ysis revealed that all these solutions are unstable with respect to linear and nonlinear spherically 
symmetric perturbations and showed that the perturbation causes the wormholes to either decay to 
a Schwarzschild black hole or undergo a rapid expansion. Here, we consider charged generalization 
of the previous models by adding to the gravitational and ghost scalar field an electromagnetic one. 
We first derive the most general static, spherically symmetric wormholes in this theory and show 
that they give rise to a four-parameter family of solutions. This family can be naturally divided into 
subcritical, critical and supercritical solutions depending on the sign of the sum of the asymptotic 
masses. Then, we analyze the linear stability of these solutions. We prove that all subcritical and 
all critical solutions possess one exponentially in time growing mode. It follows that all subcritical 
and critical wormholes are linearly unstable. In the supercritical case we provide numerical evidence 
for the existence of a similar unstable mode. 

PACS numbers: 04.25.Nx, 04.40.-b, 04.25.D- 

I. INTRODUCTION 

Wormhole spacetimes in Einstein's gravitational theory have received considerable attention in the literature. Pre- 
sumably, this is due to their interesting topological and causal properties which open the door to spectacular phe- 
nomena such as interstellar travel and time machines [H, 0, [1] ■ More recently, wormholes have also been proposed 
as black hole mimickers 3 (sec also @). However, there are several problems which make these solutions somehow 
unattractive from a physical point of view and pushes them to the verge of science fiction. First of all, wormholes need 
to be supported by exotic matter if they are to be considered as asymptotically flat, globally hyperbolic spacetime 
solutions of Einstein's field equations P, @. This means that they require an energy-stress tensor which violates 
the (averaged) null energy condition, a phenomena which has not been observed for classical matter fields in the 
laboratory. Furthermore, to our knowledge, there are no static, asymptotically flat wormhole models which have 
been shown to be linearly stable with respect to arbitrary linear fluctuations of the metric and matter fields in the 
given model. This is to be contrasted with static, asymptotically flat black holes with a regular horizon in vacuum 
or electrovacuum spacetimes, which are known to be stable with respect to linear perturbations 0, H, 0, EL E3] • 
Therefore, even if some exotic form of matter could be found in the Universe, it is not clear whether or not it could 
be used to form static wormholes. 

In previous work [l3l . ITij we analyzed the question of wormhole stability for a very simple matter model which 
consists of a massless ghost scalar field, that is, a massless scalar field whose kinetic energy has a reversed sign. 
We found that all static and spherically symmetric wormholes in this theory are unstable with respect to linear 
and nonlinear perturbations. Each of these wormholes possesses a single unstable mode which causes the wormhole 
to collapse to a black hole or to undergo a rapid expansion. Furthermore, the time scale associated to the linear 
instability is of the order of the areal radius of the wormhole's throat divided by the speed of light. 

The purpose of this article is to analyze whether or not one could stabilize these wormholes. One possible mechanism 
for stabilization is to consider stationary, rotating generalization of such wormhole solutions and hope that they become 
stable if the angular momentum is large enough. Slowly rotating wormholes have been constructed in fl~5j by studying 
linear perturbations of the static solutions. However, such slowly rotating solutions cannot be expected to be stable 
since the unperturbed solutions are unstable. Wormhole solutions of the full nonlinear field equations which represent 
rotating generalizations of the static ones have been considered in [III [l7| , but a detailed stability analysis of such 
solutions is not expected to be simple. For this reason we consider, here, a different possible stabilization mechanism 
which is the addition of an electromagnetic charge. Therefore, we add a Maxwell field to the ghost field and consider 
static, spherically symmetric wormhole solutions in this theory. It turns out that the resulting field equations can 
be integrated exactly and give rise to a four-parameter family of wormhole solutions. This is described in section HT1 
Some of these solutions have been found in |17j | . The phase space of solutions can naturally be divided into subcritical, 
critical and supercritical wormholes depending on whether the sum of their asymptotic Arnowitt-Deser-Misner (ADM) 
masses is negative, zero or positive, respectively. 

Next, we analyze the stability of such wormholes with respect to linear, spherically symmetric fluctuations. In 
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section Hill we derive the perturbation equations and cast them into a constrained wave system for two gauge-invariant 
quantities. We then discuss different ways of decoupling this system. In section [TV] we first prove that all subcritical 
and critical wormholes are linearly unstable. We do this using the same techniques as in our previous work, namely the 
theory of Schrodinger operators on the real line. Then, we focus our attention to supercritical wormholes which turn 
out to be more interesting. In this case, we are not able to reduce the perturbation equations to a single wave equation 
with regular potential. Therefore, we analyze the numerical stability of these wormholes by numerical integration 
of the constrained wave system. The results indicate that these wormholes are linearly unstable as well, with the 
perturbations exhibiting an exponential growth of the form e /3r for large values of proper time r at the throat, where 
f3 > 0. Interestingly, we find that if the charge is large enough, this exponential growth is accompanied with an 
oscillating factor of the form cos(cjt — 6) for some frequency u> and phase 8. Our results also indicate that the growth 
rate (3 decreases monotonically when the charge increases, and for the range of parameters used in our numerical 
simulations we observe that (3 can be decreased by a factor of more than 100 when compared to the uncharged case. 
Finally, we show by a numerical matching algorithm that these asymptotic solutions correspond to eigenfunctions of 
the spatial perturbation operator corresponding to a complex eigenvalue —{(3 + iu) 2 . 

A summary of our results and conclusions are given in section |V1 Technical properties of the static wormhole 
solutions which are needed for the stability analysis are stated and proved in an appendix. 



II. STATIC, SPHERICALLY SYMMETRIC CHARGED WORMHOLES 

We consider a spherically symmetric gravitational field which is coupled to a massless ghost scalar field $ and to 
an electromagnetic field F. Therefore, choosing suitable local coordinates t, x, ip, the metric has the form 

ds 2 = -e 2d dt 2 + e 2a dx 2 + e 2c (dd 2 + sin 2 dip 2 ) , (1) 

where the functions d = d(t,x), a = a(t,x) and c = c(t, x) depend only on the time coordinate t and the spatial 
coordinate x. We assume that $ and F are also spherically symmetric which means that $ = $(t,x) and that F has 
the form 

F = a dt A dx + (3 dQ A sin i? d(p, (2) 

with two functions a = a(t,x) and (3 = (3(t,x). We are interested in traversable wormhole geometries which consist 
of a throat connecting two asymptotically flat ends at x — > +oo and x ~ > -co, respectively. This means that 
the areal radius r = e c is strictly positive and proportional to \x\ for large \x\ and that the 2-manifold (M,g) = 
(K 2 , —e 2d dt 2 + e 2a dx 2 ) is regular and asymptotically flat at x — > ±oo. 
The equations of motion are 



Ruv — Ko 



F a F — - n F Gp F 



kV m $-V„$, (3) 



V"J^„ = 0, V [CT F H = 0, (4) 
VV^ = 0, (5) 

where R^ v and denote, respectively, the Ricci tensor and the covariant derivative associated with the spacetime 
metric g^ v . In terms of Newton's constant G, the coupling constants kq and k are given by k$ = 2G > and 
k = —8irG < which is negative due to the fact that <& describes a ghost scalar field. For the spherically symmetric 
ansatz (|H2p Maxwell's equations (|4]) imply that a = Q e e a+d ~ 2c and (3 = Q m with Q e and Q m two constants 
representing, respectively, the electric and magnetic charge. Setting Q := \/Q 2 + the remaining Eqs. (|3l5p yield 
the evolution equations 



d t (e a - d a t ) - d x (e d - a d x ) - e a ~ d c 2 t + e d - a c 2 x - e a+d - 2c = -K Q 2 e a+d - ic - | [e a ~ d ^ 2 - e d ~ a $ 2 ] , (6) 

d t {e a - d+2c c t ) - d x {e d - a+2c c x ) = ~e a+d + ^ Q 2 e a+d - 2c 1 (7) 

d t (e a - rf+2c $ t ) - d x (e d - a+2c $ x ) = 0, (8) 
which are subject to the constraints 

H := e d - a [2c xx + (3c x - 2a x )c x ] - e a - d c t (2a t + c t ) - e a+d - 2c + ^ Q 2 e a+d - 4c + | [e^ 2 + e^* 2 ] = 0, (9) 
M := 2ct x + 2c t c x -2d x ct-2a t c x + K$t®x = 0- (10) 
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Here, the subscript t and x refer to the derivatives with respect to t and x, respectively. 

For a static configuration, the scalar field $ and the metric coefficients d, a and c are independent of t. In this case 
the field equations can be integrated analytically. For zero charge the corresponding solutions have been obtained in 
fl8l. f]~9l| . Here, we generalize their solutions for the charged case. In the static case, the field equations simplify to 



[e d - a+2c d x ] x 
[e d - a+2c <P x ] x 



^Q 2 e a+d - 2c , 

£ a+d _ ^£q2 £ a+d-2c 

o, 



,2(o-c) _ ^£q2 e 2(a-2c) + 



(11) 

(12) 
(13) 
(14) 



Adopting a gauge where a = — d, the first two equations imply that [e 2( - c+<i ^] xa: = 2 which has the general solution 
£ 2{c+d) _ x 2 _|_ 2a\x + «o with two constants ao and a.\. By a suitable translation of the coordinate x it is always 
possible to obtain o,\ = 0. Furthermore, since we are interested in wormholc geometries with the properties described 



below equation (fl]) we need e 2(c+d) > for all x £ R. Therefore, we have e 



2(c+d) 



b with some strictly positive 



constant b > 0. Equation (TlB")) then gives $ = $i arctan(x/6) + <&o with two integration constants <&o and $x- Since 
only t he gradien t of $ appears in the equations we set the parameter $0 to zero in what follows. Next, setting 
e := ^/n Q Q 2 /2b 2 and y := arctan(x/6) £ (— k/2, tt/2), Eq. (TTT]) gives 



d = e 2 e 2d 
yy 



The unique local solution with initial conditions d\ y=Q = 70 £ M and rfj/I^Q = 71 £ K is 



d = 70 - log 



cosh(Ay) — 71 



sinh(Ay) 



A 



A 



\Jll~ e 2 T°e 2 . 



(15) 



(16) 



We distinguish between the subcritical case where A > 0, the critical case where A = and d = 70 — log(l — 
7iy) and the supercritical case where A = i\x for some real, strictly positive number in which case d = 70 — 



log cos(^y) - 71 



gin(W/) 



In order to obtain a global wormhole solution, we need the expression inside the square 



brackets to be strictly positive for all y £ [—tt/2, tt/2]. This is the case if and only if 

tanh(A 2 ) ^ ^ ^ ^ e subcritical case, 
f |7i| < 1 in the critical case, 



fj, < 1 and tan ^ 2 '' |7i| < 1 in the supercritical case 



> . 



Finally, Eq. ([14"]) yields the relation 



between the parameters A and $1. 
Summarizing, we obtain the solutions 



n$l = 2(1 + A 2 ) 



ds 2 



„2<* 



dt 2 + e" 



-2(/ 



Q m d$ A sin$ dyj, 

[dx 2 + (x 2 + b 2 ) (dd 2 + sin 2 § dip 2 )} , 



(17) 

(18) 

(19) 
(20) 
(21) 



sinh(Ay) 
'1 A 



with A = v/ 7l 2 - KQ (Q 2 + Q 2 m )e 2 io/(2b 2 ). The 



where y = arctan(x/6) and e 2d = e 270 cosh(Ay) 

parameters b, $1, Q e , Q m , 70 and 71 are subject to the two constraints (fT7|) and <|T~8|) . Notice that the constant 
rescaling i 1— > exp(— fl)t, x t~ > cxp(f2):c, 6 1— > cxp(f2)&, 70 ^ 7o + 7i — * 71) Qe l— * Oe, Q m >— > Q m with £7 a 
nonvanishing constant leaves the solution unchanged. In particular, we can rescale the coordinates such that cither 
lim d = or lim d = which shows that the spacetime described by the metric (|21l) has indeed two asymptotically 

x — > + oo x — > — 00 

flat ends at x — ► +00 and x — + — 00, respectively. 

Therefore, we obtain a four-parameter family of wormholc solutions characterized by the scale invariant quantities 
B := be~ JO > 0, 71 £ R, Q e £ K and Q m £ M, which are subject to the restrictions (TTT]) . In the particular case 
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Qe = Qm = this family reduces to the bi-parametric solution obtained in [lU, [nj. As shown in [l3|, [l4[ these 
uncharged wormholes are unstable with respect to linear and nonlinear perturbations. 

Let us analyze the physical properties of the wormhole solutions. First, the wormhole throat is given by the global 
minimum of the areal radius r = \/x 2 + b 2 e~ d . In Lemma [T] in the appendix we show that r has a unique minimum 
which is determined by the unique root of the function c x given in Eq. (|22[) below. Next, we compute the Misner-Sharp 
mass function [2pj| . For the spherically symmetric spacetime metric given by equation |T]) it is defined by 



T 6 

m (t, x) := - [1 - g(dr, dr)} = — 



1 + e 2 ^c 2 - e 2 < c 



Specialized to the static family of solutions described by Eqs. (|19I20I21[) it yields 

where c x is the derivative of the logarithm of the areal radius, 

71 cosh(Ay) 



m(x) = - [1 



b 2 



tan(y) — A 



A sinh(Ay) 



A cosh(A?/) — 71 sinh( Ay) 



(22) 



The ADM masses of the two asymptotically flat ends can be computed by considering the asymptotic values m± oc := 
lim m(x) of the mass function, which yields 

x — >±oo 



= B 



TT 



7i cosh ( A— - A sinh A— 



7T 



= —B 



71 



cosh ( A- 



A sinh A 



(23) 
(24) 



In particular, we have the relations 



m +c 



= 2B7iCosh[A- 
= -2BAsinh(A? 



m +0O m_ 



= -B 1 



' sinh" 



A^ 
2 



where we have set v := y/n (Q 2 + Q 2 n )/(2B 2 ). From the first relation and Eqs. (|19|20I21|) we see that the wormholes 
arc reflection-symmetric about their throat if and only if 71 = 0. Since A = y/'y 2 — v 2 the asymmetry parameter 71 
and the dimensionless charge v determine the asymptotic masses m+00 and to-od, the total electromagnetic charge 
Q ■= \/Q 2 + Qm an d the areal radius of the throat up to the scale factor B. From the second and third relations 
we see that in the subcritical case (A > 0) the two masses have opposite signs and that their sum is negative. In the 
critical case (A = 0) the sum of the two masses is zero, and the masses are different from zero unless 71 = Q e = Q m = 
in which case m +oc = m_oo = 0. In the supercritical case (A = ifx, \i > 0), the sum of the two masses is positive, 
with opposite signs if 7J > fi 2 tan 2 (/z7r/2) and equal signs if 7 2 < /j, 2 tan 2 (/x7r/2) while one of the masses is zero and 
the other is positive if 71 = ±/ztan(/U7r/2). 



III. DERIVATION OF THE PULSATION AND MASTER EQUATIONS 

In this section we derive the relevant equations for analyzing the linear stability of the four-parameter family of 
static wormhole solutions discussed in the previous section. For this, we consider small perturbations of the form 

$(A) = <I> + \S<f> + 0(\ 2 ), 

where $ is the background solution, and where 

5$ := A$( A ) 

denotes the variation of (f>. The same applies to the other fields d, a and c. A general method for analyzing such 
perturbations has been developed in (2lj | . Since in spherical symmetry there are no gravitational nor electromagnetic 
dynamical degrees of freedom, one obtains a single master equation for the linearized scalar field <5$. However, as we 
have discussed in detail in [l 31 ] for the uncharged case, the resulting master equation turns out to be singular at the 
throat which leads to difficulties when studying the stability by standard methods based on Schrodingcr operators. 
For this reason, we will not use the method described in [2l| and instead base our treatment on a gauge-invariant 
approach which leads to a constrained wave system which is everywhere regular. 
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A. Gauge-invariant quantities 

With respect to an infinitesimal coordinate transformation St i— ► 6t + £*, <5x i— > <5x + £ x on the 2-manifold M 
generated by a vector field (£ t ,£ x ), we have 

<Sa^<5a + e - a (e a f% , Sc^Sc + Cc x , 6$ t-> 6$ + £ x $ x . (25) 
Since $> x ^ everywhere we may construct the following two gauge-invariant fields, 

A := Sa - e~ a [ e a — , C := Sc - c x — , (26) 

which reduce to Sa and Sc, respectively, in the gauge <5<& = 0. 

B. Constrained wave system for A and C 

Here we derive a constrained wave system for the gauge-invariant quantities A and C. In order to do so, we first 
rescale the coordinate x such that 6=1. In terms of the coordinate y = arctan(a;) which satisfies d y = e d ~ a+2c d x , 
the background Eqs. (TTlJfjJJ) yield 

d yy = e*e™, c yy =e 2 ^~e 2 e 2d , <f yy = 0, (27) 

where e 2 = kqQ 2 /2. 

Next, we consider the constraint equations (|9I10|) . Their linearization yields 

5(e a - d H) = 25c xx + {6c x - 2a x )5c x ~ 2c x 5a x - 2(Sa - Sc)e 2(a ~ c ^ + n Q 2 (5a - 25c)e 2a - 4c + K<f> x 5$ x = 0, (28) 
SM = 2Sc tx + 2{c x ~ d x )8c t - 2c x Sa t + k $ x <5$ t = 0. (29) 



With the help of the background equations (|27|) one may rewrite this as 
1 



^e 2{d ~ a+2c) 8{e a - d -H) = 5cy + (c v -d y )5c-c y 5a+~$ y 6$ =0, (30) 



2 



i e d-a+2c SM = 

2 



5Cy + (Cy - dy)SC - CySCL + ~<f>yS$ =0, (31) 

2 J t 



which shows that the expression inside the square bracket must be equal to a constant a. In terms of the gauge- 
invariant quantities A and C defined above, the resulting first integral is 

Cy + (Cy - dy)C ~ Cy A = G. (32) 

The interpretation of the constant a is the following. With respect to an infinitesimal variation of the constants -B71 
and A (keeping the charges Q e and Q m fixed), the family of static solutions (|2Tj) yields the linearized solution 

C=^ 2 8{B 11 )-F 5 -±, A = -(l + xy)^^+C, (33) 



where the function F is given by 



f = 1 + IS _^ = 1 + <^. (34) 



Introducing the expressions (f33|) in ([32]) gives a = —S(Bji). Therefore, a describes variations of the static family of 
wormhole solutions with respect to the constant B71. Since any solution to the linearized equations may be written 
as the sum of such a variation plus a solution with 8{B^ii) = we may assume that a = in the following. 

Next, we linearize the evolution Eqs. (|6I7I8[) . For simplicity, we choose the gauge such that 8$ = 0, in which case 
A = Sa and C = 8c. Linearization of Eq. §8$ yields 8d — 8a + 28c = h(t) for some function h(t) which we may set 
to zero by a redefinition of St. Using 8d = Sa — 25c and the first integral (|32|) with a = in the linearization of the 
evolution Eqs. (|6|7|) a lengthy calculation yields the constrained wave system 

u tt - e~ 2c [e~ 2c u y ] + Vu = 0, C := (u 2 ) y + (c y - d y )u 2 ~ c y u x = 0, (35) 
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where we have defined 

A - G \ 4c / 3c 2 + Ac y d y - 3e 2 ^) + 5e 2 e 2d 4A 2 \ 

CJ> V ~ 6 { -Ac 2 + 2e 2 ( d+c ) - 2e 2 e 2d 3c 2 - Ac y d y - e 2 ( d+c > + 3e 2 e 2rf ) ' [6b) 

The constrained wave system (|35j) describes the dynamics of the two gauge-invariant linearized fields A and C. The 
linear stability properties of the wormholes are determined by the Cauchy evolution of this system. Two difficulties 
with analyzing the properties of the solutions are the fact that we are confronted with a coupled system of two 
equations (as opposed to a single, scalar equation) and the presence of the constraint C = 0. In the next subsection 
we start by deriving a decoupled equation for the constraint field C based on a factorization method. This shows that 
it is sufficient to enforce the constraint and its time-derivative at an initial time. As a byproduct of our factorization 
method, we also obtain a master equation for a quantity V\ defined below, from which U\ and u 2 can be reconstructed. 
However, this equation turns out to be singular at the throat, and as mentioned before, this means one has to be 
careful with the stability analysis. For this reason, we derive in the following subsection a different master equation 
which, in the critical and subcritical cases is everywhere regular and allows to prove that such wormholes are linearly 
unstable. In the supercritical case, however, both master equations turn out to be singular, and so the constrained 
wave system (|35p has to be analyzed directly. If the potential V would be symmetric, or if it could be brought into 
symmetric form by a linear transformation of u, one could analyze the system by spectral analysis of the formally 
sclf-adjoint operator 7i = —e c d v e c d y + V. The transformation w% = U\, Wi = —U\ + u 2 brings the system into 
the form 

Wtt - e"- [e-^ Wy ] y + V W = 0, V = e~*° ( D + \ D % 

with D = 3c 2 , - 2e 2 ( c+d ' +4e 2 e 2d and E = Ac y d y - e 2( - c+d '> + e 2 e 2d + 4A 2 . If A 2 > this can be symmetrized by a trivial 
rescaling of w; for A 2 < 0, however, V12 and V21 have different signs. In fact, we will show in the next section that in 
the latter case the operator H may have complex eigenvalues, and so it cannot be written as a symmetric operator. 

C. Factorization of the Hamilton operator and master equation I 

Consider the two-channel Schrodinger operator 

H := -d 2 +V, d:= e - 2c d y = e d ~ a d x . 

Here, we try to factorizc it in the form TL = AB, with the two first-order operators 

A = d + K, B = -d + K, 

where K is a 2 x 2 matrix which has to satisfy the Riccati matrix equation 

dK + K 2 = V. (37) 

If K solves (|37|). the factorization TL = AB allows us to rewrite the wave problem uu + TCu = into first-order form, 

u t = Av, (38) 
v t = -Bu. (39) 

In particular, it follows that v satisfies the dual wave problem 

v tt + BAv = 0, (40) 

where BA = —d 2 + V with the transformed potential V = —V + 2K 2 . 

A solution to (|37|) can be obtained as follows. We demand that the second component of v t = —Bu is proportional 
to the constraint variable C defined in Eq. ([33)1 . This implies that K must be of the form 



= e " C ( I d -C, ) 
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with two unkown functions / and g. (Notice that in this case dtV2 = e 2c C.) Introducing this ansatz into Eq. (f37|) 
yields the unique solution 



f = dy 

for / and g. The transformed potential is then 



1 + A 2 




2(1+A 2 ) 



A 2 

.9 = 



_ _2 _ 2A- 



(42) 



Therefore, the constraint variable i>2 satisfies a decoupled wave equation with potential V22 = e (c vy — cj), and it 
is consistent to enforce the constraint C = 0. Setting «2 = the first-order system (|38I39[) then reduces to 



dtUi = dv\ + e~ 2c fvi , 

d t U2 = e~ 2c c y vi , 

d t vi = dux - e~ 2c (fui + gu 2 )- 



In particular, v\ satisfies the following master equation, 



d 2 -d 2 + V u 



vi = 0, 



Vu 



-4c 



2(1 + A 2 ) 



(43) 
(44) 
(45) 



(46) 



As mentioned above, the resulting potential Vu is singular at the wormhole throat, where c y = 0. As discussed in 
detail in [l3[ this enforces an unphysical boundary condition at the throat if one tries to define — d 2 + Vu as a self- 
adjoint operator. Namely, it requires v\ to approach zero sufficiently rapidly as y converges to the throat's location. 
However, there is no reason for enforcing such a strong condition on v\ from a physical point of view. As we will sec 
in the next section, physically permissible perturbations even allow v\ to diverge at the throat. In particular, this 
implies that the operator —d + Vu is not symmetric when defined on the space of physically permissible states. 

As we show next, it is possible to obtain a different master equation which in the subcritical and critical cases is 
everywhere regular and yields a self-adjoint operator without enforcing unphysical boundary conditions. 



D. Master equation II 

The derivation of the new master equation is based on the observation that the constrained wave system ([35 
possesses the particular solution 



^static 



G \-,-cf 1^(1 + 1 (47) 



H ~ \ F 



which is obtained from ((33)) after setting S(B-fi) = and 5 A = —A. A related solution is obtained by multiplication 
with t, corresponding to the following solution of the first-order system (|43I44I45[) 

U = tu StatiC , Vl = * = — ■ 

Cy 

Since is a time-independent solution of the master equation (|4"6")l we may rewrite the latter in first-order form 

d tVl = (d+^jx, (48) 

The corresponding dual wave equation is 

[d 2 - d 2 + W] X = 0, (50) 
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with the transformed potential 

_ e _4c J 0/1 i A 2 



G r , .9-1 n 9 /G^ 



-3(1 + A 2 ) + 2d w (d„ - c y ) - 4- [c,rf y + 1 + A 2 ] + 2c 2 (^-J j , (51) 
where we have used the identity d^ 1 " + / = g c?/ . I n terms of the new variable x the first-order system yields 



dtui = d tX + e- 2c —(cyVi), (52) 
rl 



d t u 2 = e- 2c (c yVl ), (53) 

d t (c y v{) = c y d X + e- 2c (^c 2 y -c y d y -l-A 2 ^x, (54) 

which allows to obtain the gauge- invariant perturbation quantities U\, u 2 (and c y V\) from x after a time integration. 

In contrast to the first master equation, the new master equation (|50[) is regular at the throat. In fact, the potential 
W is everywhere regular as long as the function F does not have any zeroes. This turns out to be the case for the 
critical and subcritical cases, see Lemma [2] in the appendix. 



IV. LINEAR STABILITY ANALYSIS 



In this section we discuss the linear stability of the wormhole solutions in the subcritical, critical and supercritical 
cases. In the first two cases, we show instability by proving that the master equation ([50)) which is regular in those 
cases possesses precisely one exponentially in time growing mode. In the supercritical case, both master equations 
are singular, and we analyze the stability by different means. 



A. The subcritical case 



The results from the previous section allow us to describe the stability problem in the subcritical case by the regular 
master equation ([50| on the Hilbert space X = L 2 (R, e~ 2d dx) which admits the zero mode 

Since the function c y /e c is uniformly bounded on —00 < x < +00 and possesses exactly one zero, and since the function 
F is strictly positive and satisfies F/x — > A 2 (l + A 2 ) _1 7r/2 for x — » ±00, this mode belongs to X and represents a 
bound state of the Schrbdinger operator — d 2 + W. Because it has one node, it follows from the nodal theorem 1 that 
it is the first excited bound state and so the operator — d 2 + W possesses precisely one negative eigenvalue — (3 2 < 
with eigenfunction xp, corresponding to an exponentially growing mode of (|50|) which is of the form 

X{t,x)=ef }t x (x). (56) 

Since the coefficients in the Eqs. (|52I53|54|) are everywhere regular this gives rise to a unique unstable mode for each 
subcritical wormhole. We conclude that all such wormholes arc unstable with respect to linear, spherically symmetric 
perturbations. 



B. The critical case 

Linear perturbations in the critical case arc also described by the regular master equation (|50p . However, in contrast 
to the previous case, the function 

x{l - iw) - 7i 



Xo = l/*o 



Vl + x 2 {l - 7iy) 2 



1 See, for instance, or [23ll for a generalization to systems. 
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is not normalizablc, so the above argument based on the nodal theorem docs not directly apply. Instead, we construct 
a family /„, n = 1,2, 3, of static solutions to (|50[) which are defined on the interval — n < x < +oo and satisfy the 
following two properties for large enough n: (i) f n (—n) = 0, (ii) (f n )x(—n) ^ 0, (iii) /„ has exactly one zero on the 
interval x > n. It then follows from the results in Ref. [23| that there is a unique bound state with negative energy, 
as in the subcritical case. Therefore, all critical wormholcs are linearly unstable as well. 

The family /„ is defined as follows. Let xthroat be the value of x at the throat, where c y = 0, and let n > —xthroat- 



Then, 



fn(x) := xo(x) 




, 71 ^ X <C X ' throat i 
3 ^> Xthroat ; 



where k is a constant to be determined. It is simple to verify that f n satisfies the relation Xodf 7 i = 1 + fndxo an d the 
master equation (|50[) on the two open intervals (— n, xthroat)-, (%throat, +00), and that f n (—n) = and (f n )x{~ n ) < 0. 
In order to analyze the behavior of f n near x = xthroat we first introduce the new coordinate 

p = R{x):= / e~ 2<i(s) c&, -oo < x < +00, (57) 



in terms of which we have d = d p . Next, since xo satisfies (— d 2 + W)xo = and vanishes at p = 0, we can write it in 
the form xo(p) = a.p[\ + p 2 q{p)] 1 where q is a smooth function on R such that q(0) ^ 0, and a > 0. Also, 1 + p 2 q(p) 
has to be strictly positive since otherwise xo would have more than one zero. In terms of this, we find 



Up) 



a 



R(n) _ 
P n f 2g(p)+p g(p) j- ^27 ^ n 



which shows that f n is continuous at p = 0. Furthermore, we see that (/ n ) P is also continuous at p = provided we 
choose k such that 

R(n) 

1 1 f 2q(p) + p 2 q(p) 2 



01 k R(-n) R(n) J [l + p 2 q(p)} 2 d(> ' 

R(-n) 

Therefore, /„ can be extended on the whole interval — n < x < +00. Finally, we notice that /„ is negative on 
the interval — n < x < x t h roat and positive for large enough x > x t hroat- Hence, /„ has at least one zero at some 
x* > Xthroat- On the other hand, because xo(x*)df n (x*) = 1 + f n (x*)dxo{x*) = 1 and xo{x*) > it follows 
that df n (x*) > which means that this zero is unique. In the reflection-symmetric case where 71 = we have 
Xo{x) = x/y/1 + x 2 and f n (x) = (x + n)(x - l/n)/\/I + x 2 . 



C. The supercritical case: Numerical integration of the constraint wave problem 

For the supercritical case, F has always two zeroes and both master equations are singular. Therefore, the previous 
arguments cannot be used to establish the linear instability of the wormholcs in this case. Short of an analytic proof, 
we shall analyze the stability of supercritical wormholcs by numerical means. We start in this subsection with a 
numerical integration scheme of the constraint wave problem (|35|1 . We model this scheme on the recent work in Ref. 
[23 | . where hyperboloidal time slices with a compactified space coordinate are used. One of the main advantages of 
this scheme is that the Cauchy evolution is performed on a compactified domain where the boundaries of the domain 
correspond to future null infinity. In this way, one avoids the problem of introducing an artificial timelikc boundary 
with absorbing boundary conditions. This is particularly attractive for our stability problem since we do not want 
the time evolution of our solution to be contaminated by artificial boundary conditions. 

The starting point of the numerical scheme is to rewrite the wave equation in Eq. (|35[) in geometric form 



- g ab V a V b u + Vu = 0, 



(58) 
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where g = —e 2d dt 2 + e~ 2d dx 2 is the wormhole 2-metric, V is the associated connection, and V — e~ 2d V is the rescaled 
potential. An observation of later importance is the fact that the rescaled potential decays at least as 0(x~ 2 ) when 
x — > ±00. 

Next, we introduce a new time coordinate r := t — h{x), where the height function h, which will be specified later, 
is such that the r = const slices are everywhere spacelike and asymptote to outgoing null geodesies as x — > ±00. 
This means that the function e 4d h 2 is strictly less than one and converges to one as x — * ±00. The r = const slices, 
as embedded in the four-dimensional, spherically symmetric spacetime with 2-metric g and areal radius r have mean 
curvature 

c = ~^a~ ' (59) 



where the function J is given by J := c zd h x j ^/\ — e 4d h 2 . Here, we choose the height function such that J = 
cqx(x 2 + 3b 2 ) /(x 2 + b 2 ) with some positive constant Co- This implies that the mean curvature of the r = const slices 

is c = Co 1 — 2 4f aL x x -i+b i which converges to cq > as x — ► ±00, and the function e 2d h x also satisfies the properties 

described above Eq. |59|) . 

In a next step, a new, compactified space coordinate z G [—1,1] is introduced which is related to the coordinate 
x G (—00, +00) via the transformation x = z/VL(z) with f2 € C°°[— 1, 1] a smooth function which is strictly positive 
on (—1, 1), vanishes at the endpoints z = ±1, and which satisfies the inequality L := — ztt z > everywhere on 



T, 1]. Here, we choose Q(z) := 1 



z 



2 



With respect to the new coordinates (t, z) the 2-metric assumes the form 

g = fT 2 <7, g = -a 2 d,T 2 + ^ 2 (dz + fidr) 2 , 

where 



n 2 e 2d + (nj) 2 , 7 = 4, = 

a 7 

Since f2J is everywhere regular and positive near z = ±1, the conformal 2-metric g is regular for all 2 e [—1, 1]. Since 
the two-dimensional wave operator is conformally covariant, the wave equation f|58[) is equivalent to 



- g ab V a V b u + Vu = 0, (60) 

where V = Q,~ 2 V = (x/z) 2 V is everywhere regular on z E [—1,1]. For the numerical implementation, we cast ([60]) 
into first-order symmetric hyperbolic form for the six fields u = (u±, 112), D = (-Di, -D2), n = (IT, II2), 

u T = ciTI + 7/iD, (61) 

D T = l(&n + ^(3D) z , (62) 
7 

n T = i((xD + 7/3n) z -aVu. (63) 
7 

This system is to be integrated on the compact interval z 6 [—1, 1]. The characteristic speeds arc 



d a 1 
A± = ±-+0= - 

7 7 



±jQ 2 e 2d + (nj) 2 - nj 



(64) 



Therefore, A_ < < A_|_ on z £ (—1, 1). At the left boundary, A_ = while at the right boundary, A + = which 
implies that both boundaries are outflow. This is of course expected from the fact that the boundaries z = ±1 
represent future null infinity. 

In terms of the compactified coordinates, the constraint C = reads 

Q 2 C = (Q 2 + z 2 ) e - 2d (aD 2 + 7/3II 2 ) + n 2 (c y - d v )u 2 - n 2 c yUl = 0. (65) 

A simple way of solving the constraints is to specify initial data for Tl 2 and D 2 which are compactly supported 
away from the throat, and to use the constraint C = and its time derivative in order to determine the remaining 
fields iti, 111 and D\, keeping in mind that 

III = 4 (d T Ui - 7/3Di) , D 1 = \d z u x . 
a V / 7 

The procedure we used to study the solution is as follows. 
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1. We set up initial data by specifying an initial pulse for U2 and solving the constraint, as described above. 

2. We evolve the perturbation u using the system of equations (|61I62I63[) . 

3. We measure the amplitude of m at the location of the throat. 

4. We fit the resulting value of u\ at the throat using the ansatz Acos(ujt — <5)e' 9 ' r , where r is proper time at the 
throat, and determine in this way the growth rate (3 of the perturbation and the frequency of oscillation cu in 
case there is one. 

We use a finite differences approximation method with a method of lines for the evolution of the perturbation 
(|61l62l63p . When the dimensionless charge v is not too close to its limit value given by the constraint (|17p . second- 
order accurate stencils show convergence and confident results. In the case v approaches its limit value, the results 
obtained from the second-order accurate stencils fail to show convergence, and we use eighth-order accurate stencils 
in order to reduce the errors and work in the convergence regime. 

In Fig. 1 we present our results for the reflection-symmetric case 71 = 0. For small values of v the perturbation 
grows exponentially and no oscillating mode shows up. However, there is a threshold of v around v* = 0.55 above 
which the exponential growth is modulated by a harmonic component. Also, we observe that the growth rate decreases 
to less than 1% of the rate in the uncharged case as v increases from zero to 0.96. Above v = 0.96 the functions in 
the evolution equations become stiff and our numerical approach breaks down. In order to validate our results, we 
compare them with the results of the matching method described below. 

An interesting question that arises from the plots in Fig. 1 is whether or not the growth rate (3 may be zero for 
some value of v lying between 0.965 and 1, implying the existence of charged wormholes which arc stable with respect 
to linear radial perturbations. A different possibility is that /3 stays positive for all < v < 1 and converges to zero 
for v — ► 1, meaning that all supercritical wormholes are linearly unstable, but the growth rate can be made arbitrarily 
small by adding a sufficient amount of charge. The third possibility is that (3 is bounded away from zero for all 
< v < 1 in which case all supercritical wormholes are linearly unstable as well. The answer to this question requires 
an analytic understanding of the behavior of f3 as a function of v near one and lies beyond the scope of this paper. 



S 0.6 



v=0.96 



Evolution method 2nd order 
Evolution method 8th order 
Matching method 



small values of v 



Evolution method 8th order 
Matching method 



0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 
V 



FIG. 1: Left panel: values of the frequency ui and growth rate /3 of the perturbation for the reflection-symmetric case 71 = 
and several values for the dimensionless charge v using three different techniques. Each point in this plot corresponds to a 
given value of v. The results in this plot indicate that the exponential growth rate decreases as the charge v is increased. Right 
panel: the growth rate fi versus v for large values of v. 



In Fig. 2 we present results for the asymmetric case for different values of 71 7^ 0. The behavior we find is similar 
to the massless case: for each 71 there is a threshold value for v above which the perturbation shows an oscillatory 
behavior while growing exponentially. 



D. The supercritical case: Eigenvalues of the pulsation operator 



The numerical results from the Cauchy evolution suggest the existence of eigenmodes of the pulsation operator 
with time-dependency e (l 3 + w ) t . If suc h a n eigenmode exists, it must satisfy the constrained wave system (|3"5j) . with 
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FIG. 2: Values of the frequency and growth rate of the perturbation for different values for 71 and v. In this case the allowed 
values of v are restricted by the condition (I17|l for the supercritical case. Similarly to the reflection-symmetric case there is a 
threshold between purely exponential growth and exponential oscillating growth. 



the functions ui(t, •) and u 2 (t, •) being regular for all times t > 0. Defining the quantity "J" := e c c y v\ = e c dtU2 and 
using master equation I, we may rewrite the system (|43l44l45p in the form 



d t ui 
d t u 2 



1 



2e c c 



mi 



[* m - 2d yv * - e 4c *„] , 



Cy{U\)y — (Cydy + 1 + A )U\ + A^] 



(66) 

(67) 
(68) 



where we notice the fact that c yy is everywhere positive, see Lemma [T] in the appendix. Since all coefficients in these 
equations are regular, it follows that ^> is regular if and only if u\ and u 2 arc. Therefore, we look for eigenmodes for 
which "J" is regular and has time-dependency eS ,3+luJ ' >t . 

As a consequence of master equation I, the quantity "J" satisfies the wave equation 



e 4c * t 



2c„ 



2 d 



* = 0. 



(69) 



We are looking for solutions of this equation with time-dependency e^ + ' iw ) t which are everywhere regular and vanish 
in the asymptotic regime y <— ► ±7r/2. This leads to the following eigenvalue problem, 

d,, 



- * 



2c 



uu 



!l!l 



llll 



# = -e 4 T 2 *, 



(70) 



where T = /3 + iuj. This equation has a regular singular point at the throat c y = 0. In the reflection-symmetric case 
where 71 = one finds 



2c ? i 



l + -(l + V 2 )y 2 + 0(y 4 ) 



L yy 



4 2 2 



0(y 4 ), 



Ac 



1 + 2(1- v 2 )y 2 + 0(y 4 ), 



and wc obtain a reflection-symmetric, local solution of the form 



1. 



2„,2 



*l(l/) = l- o r v 



1 

24 



r 3 p4 + ^2(5^2 _ j j + ^2] y 4 + 0(y6) _ 



(71) 
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In order to understand the asymptotic behavior of the solutions of Eq. ((70)) when y — > ±ir/2 approaches the asymptotic 
regime, it is easier to go back to the consideration of master equation I, which in terms of the coordinate p defined in 
(1571) reads 



d;vi + 



Vn(p) 



vi = 0, 



(72) 



where we have assumed again a time-dependency e for v\ . In the reflection-symmetric case, we find 



2/itan(/i|) 1 

COS 4 X 3 



[ - 

X 



According to standard theorems 2 , there exist local solutions near p = +oo such that 



[i + Hp)] , 



where h is a C 2 -function satisfying h(p) — > 0, h p (p) — > when p — > +oo. Therefore, we obtain the asymptotic solution 



*2 = e- c c v e- T P [1 + h(p)} 



(73) 



for y — > 7r/2. Our numerical method consists in integrating the local solutions (|71I73|) and to match their Wronski 
determinant 

W(v)-det( * 2(y) 

at some intermediate point < y\ < tt/2 by fine-tuning the complex eigenvalue T. The result is shown in Fig. 1, and 
agrees very well (within 5% of relative error) with the growth rate and oscillation frequency found from the Cauchy 
evolution. 



V. CONCLUSIONS 



In this article we analyze the stability of static, spherically symmetric general rclativistic wormhole solutions sourccd 
by a massless ghost scalar field and an electromagnetic field. To this purpose we first construct the complete spectrum 
of such solutions. Among the solutions we distinguish between three types depending on the values of the charge and 
mass parameters: subcritical, critical and supercritical. We show that in the first two cases all solutions are unstable 
with respect to linear spherically symmetric perturbations. This is done by reducing the perturbation equation to a 
Sturm-Liouville problem and using standard tools in Schrodinger operator theory. 

In the supercritical case we are not able to reduce the perturbation equations to a scalar equation which is everywhere 
regular. We instead obtain a constrained wave system for two gauge-invariant quantities. We study this system 
numerically as a Cauchy evolution problem based on a domain compactification scheme and observe the growth of the 
perturbations at the wormhole throat. The analysis reveals: i) for small values of the dimensionlcss charge parameter 
v the growth of the perturbation is exponential in proper time, ii) there is a threshold v* in the value of the charge 
above which the exponential growth is modulated by an oscillatory component, iii) the growth rate decreases when the 
value of the charge increases, iv) as the value of the charge approaches its limiting value the functions in the evolution 
equations become stiff and our numerical approach breaks down. However, we find consistent results up to values of 
v = 0.96. Within this regime we can decrease the growth rate by a factor of more than 100 compared to the uncharged 
case. To validate our numerical results we also analyzed the eigenvalues of the spatial operator in the perturbation 
equation. Based on a matching method we find a real eigenvalue for v < v* . Above that value the eigenvalue acquires 
a nontrivial imaginary part, explaining the oscillatory behavior found in the Cauchy evolution. The eigenvalues fit 
the growth rate and the oscillation frequency obtained with the Cauchy evolution within an accuracy of 5%. The fact 
that we find complex eigenvalues provides an explanation for the inability to reduce the perturbation equations to a 
Sturm-Liouville problem, for which all eigenvalues are necessarily real. 

The motivation for this work was to study a possible mechanism for stabilizing the static, spherically symmetric 
wormholes solutions supported by a ghost scalar field. We have analyzed the stability behavior when such solutions 
are charged up by an electromagnetic field. The fact that this mechanism does not seem to be able to stabilize the 
wormholes raises some doubts about the success of other similar mechanisms, like adding angular momentum. 



2 See, for instance, Ref. [25t . 
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APPENDIX A: SOME TECHNICAL PROPERTIES OF THE STATIC WORMHOLE SOLUTIONS 

The family of charged, static wormhole metrics is determined by the function 



7o 


- log 


cosh(Ay) 


7o 


- log 


1 - ny] 


7o 


- log 


cos(fiy) - 



sinh(Ay) 



,A>0, 
,A = 0, 

, A = i/i,0 < fj, < 1, 



where A := \J — e 2 ~t°e 2 and y = arctan(x/6). Here, 71 is subject to the inequalities (|17|) . For simplicity, we choose 
6=1 and 70 = in what follows. Then, the metric functions d and c = log Vx 2 + 1 — d satisfy the following relations, 

d yy = e 2 e 2d , c yy = l+x 2 ~e 2 e 2d , (2dy + Cy)c y = x 2 -e 2 e 2d -A 2 , (Al) 

which follow from Eqs. pill2ll4|) with e 2(c+d ^ = 1 + x 2 , e 2 = n a Q 2 /2 and the relation dTSJ. This together with 
c y = x — d y also implies the equations 

dl=A 2 + e 2 e 2d , dyy = d y — A 2 , c yy = 1 + A 2 + x 2 - d\ (A2) 

which turn out to be useful. The next result is related to the global behavior of the areal radius r = e c . 

Lemma 1 The function c has a unique local minimum and c yy is strictly positive on the interval —00 < x < +00. 

Proof. First, we notice from Eq. (|A2jl that at points where c y = we must have c yy = 1 + A 2 > which shows that 
local extrema of c are necessarily local minima. Since c — > +00 when x <— > ±00 it follows that c has a unique local 
minimum at some point x = x throat ■ 

Next, we prove that c yy is everywhere positive. For this, we first notice that in the uncharged case e = 0, it follows 
from Eq. (|A1[) that c yy = 1 + x 2 > 1. Next, assume that e > and that c yy has a zero at some point x*. Since 
c.yy = 1 + A 2 > at Xthroat, x* lies either to the left or to the right of Xthroat- Suppose that x* < Xthroat- Then, from 
Eq. (|A2j) . it follows that c yyy = 2d yy Cy < at x = x*, where we have used d yy = e 2 e 2d > and c y < for points to 
the left of the throat. This means that c yy crosses the x-axis from above at all its zeroes lying to the left of Xthroat- 
However, since c yy > at Xthroat-, this means that c yy cannot have zeroes to the left of x t hroat- A similar argument 
shows that c yy cannot have zeroes to the right of x t hroat either. rj 

Next, we analyze the properties of the function 

F = l + xy- yCy = l+ (^ x + d y)y 
+ V 1 + A 2 + 1 + A 2 

defined in p4]) . Using Eqs. (|A2[) it is not difficult to show that F satisfies the relation 

F y = d y F + (1 + xy)-A—Cy . (A3) 
For the linear stability analysis, it is important to know whether or not F has zeroes. 

Lemma 2 In the subcritical and critical cases, the function F is strictly positive on the interval —00 < x < +00. In 
the supercritical case, F has precisely two zeroes. 

Proof. In the critical case, F = (1 — 71?/)" 1 which is everywhere positive. In the other cases, since d y is bounded and 
A 2 7^ 0, we see from its definition that F satisfies F — > +00 (F — > —00) as \x\ — ► +00 in the subcritical (supercritical) 
case. On the other hand, at Xthroat where c y = 0, we have F = 1 + xy > 0. Finally, let x* be a zero of F which lies 
to the left of Xthroat- From Eq. (|A3[) it follows that (1 + A 2 )F y = (1 + xy)A 2 c y at x* which is negative (positive) in 
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the subcritical (supercritical) case. This means that all zeroes of F lying to the left of x t hroat cross the £-axis from 
above (below). A similar argument shows that all zeroes of F lying to the right of x t hroat cross the x-axis from below 
(above). As a consequence, F has no zeroes in the subcritical case and two zeroes (one to the left and one to the right 
of Xthroat) in the supercritical case. q 
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